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Effective spin-orbit coupling can be created in cold atom systems using atom-light interaction. We 
study the BECs in an optical lattice using the Gross-Pitaevskii equation with spin-orbit coupling. 
Bloch states for the linear equation are numerically obtained, and compared with stationary solutions 
to the Gross-Pitaevskii equation with nonlinear terms. Various vortex lattice states are found when 
the spin-orbit coupling is strong. 

PACS numbers: 03.75.-b, 03.75.Mn, 05.30.Jp, 67.85. Hj 

Recently, Bose-Einstein condensates (BECs) with effective spin-orbit coupling were created in cold atom systems 
using atom-light interaction The spin-orbit-coupled BECs are actively studied theoretically Q. Wang et al. found 
that the mean-field ground state has two different phases: plane-wave and stripe phases depending on the nonlinear 
interactions Q • Half vortex states were found in a spin-orbit coupled BECs confined in a harmonic potential [3, [B| . 
Exotic spin textures were predicted in Bose-Hubbard models corresponding to spin-orbit coupled BECs in the Mott- 
insulator phase [f| 0] • 

The GRoss-Pitaevskii (GP) equation is a mean-field approximation for the BECs with the spin-orbit coupling. 
There are some studies for the GP equation with spin-orbit coupling in optical lattices [H, Q. In this paper, we study 
vortex lattice solutions to the GP equation in a square type optical lattice. The model equation is expressed as 

= -iv 2 ^ + + (5l^+| 2 +7lV'-| 2 )V' + -e{cos(2^) + cos(27ry)}^ + + A^''' 



— i- 



dt 2 ' " ' ' \ dx dy 



■ ()l -V>_ + (^-r + 7|V'+r)V'--e{cos(27rx) + cos(27ry)}^_ + A^-^-i-^J, (1) 

where if) = denotes the wave function of the spinor BECs, e is the strength of the optical lattice, g and 7 

express the strengths of interactions respectively between the same and the different kinds of atoms, and A denotes 
the strength of the Rashba spin-orbit coupling. We have assumed that the wavelength of the optical lattice is 1. 

If g and 7 are zero, Eq. (1) becomes linear equations with spatially-periodic potential. The Bloch states are 
stationary solutions to the linear equations, which are expressed as 

ip + (x, y, t) = 4>+{x, y) exp(ik x x + ik y y - i[it), ip-(x, y,t) = 4>~(x,y) exp(ik x x + ik y y - ifit), (2) 

where 0+ and <f>- are periodic functions of wavelength 1. Therefore, <j) + and 0_ satisfy 
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The eigenvalue fi and the eigen function <f> + and </>_ can be numerically obtained from the stationary solution of the 
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FIG. 1: Eigenvalue fi vs. k x for (a) A = n/2, k y — 0, (b) A = 3n/2, k y = 0, and (c) A = tt, k y = k x . Dashed curve in Fig.l(a) is 
plotted using Eq. (8) and the dashed curve in Fig. 1(c) is obtained using Eq. (6). 



linear equation [Toj : 



1 V 2 0+ -\{k 2 x + kl)(j) + + ik x ^- + ik y ^p- + e{cos(27ra) + cos(27ry)}0 + 



dt - 2 V ^ 2 \^»ViV+^>»* ^ ^ Qy 

' d<j>- .£?</>_ 
dx dy 



-A ( — — + ik x <fr- + k y 4>- ) + (!</>+, 



= ^V 2 0_ - hk 2 x + k 2 )(t>- + ik x ^— + ik y ^— + e{cos(27nr) + cos(27ry)}0_ 



dt ~ 2 V v ~ 2 yn ^^ JV -^" nx Ox ^ dy 



dfx 



a(N - TV), (4) 



where a > is a parameter and fixed to be 5 in our numerical simulation. N = J Q J^d^+l 2 + \(t>~\ 2 )dxdy is the total 
norm, and No is fixed to be 1 by the normalization condition. The time evolution of the dissipative equation (4) 
leads to a stationary state and the total norm N approaches Nq — 1. The eigenvalue [i in Eq. (3) is obtained as /1 in 
Eq. (4) at the stationary state. In this numerical method, the ground state for fixed values of k x and k y is obtained 
at the stationary state, starting from most initial conditions, because the total energy decreases in the time evolution 
of Eq. (4). Excited states are obtained by removing the ground state by the method of orthogonalization. Figure 
1(a) shows fi(k x ) as a function of k x for k y = 0, e = 5, and A = n/2. fi(k x ) is a periodic function of k x with period 
2it. There are peaks near k x = 0,7r and 2-7T and minima at k x ~ n/2 and 37r/2. The peak point at k x = 7r is a cusp 
point, where two fi(k x ) curves corresponding to the ground state and the excited state cross, although the branch 
of the excited state is not shown. For A = 0, n(k x ) increases monotonously as k x increases from and reaches the 
maximum at the edge of the Brillouin zone at k x = n. If there is no optical lattice, i.e., e = 0, fi(k) takes a minimum 

at k — A where k = \Jk x + ky [H, The minimum point for A = n/2 locates near k x = A, and the peak corresponds 

to the edge of the Brillouin zone. Figure 1(b) shows fi(k x ) at A = (3/2)tt. There are a large peak at k x = and 
2tt and a small peak at k x = tt and minima at k x = 2 and k x = 4.3. Figure 1(c) shows /i(fc) as a function of k x for 
k x — k y , A = 7T, e = 5. There is a large peak at k x = and 2ir, a small peak at k x — tt, and minima at k x ~ 2.5 and 
3.8. The wavenumber k x ~ 2.5 is close to X/y/2 <~ 2.22 by the simplest approximation k — A but slightly deviated. 
The approximation fc m j n = A for the minimum point of fJ,(k) becomes worse for large A. The small peaks correspond 
to the edge of the Brillouin zone. 

Figure 2(a) shows \(f>+(x,y)\ and \<fi-(x,y)\ as a function of y in the section x = at k x — 3tt/2 for A = 37r/2. The 
modulus |0+| and \4>—\ take maximum at different positions. The minimum value is almost zero, which implies the 
existence of vortices. Figure 2(b) shows a contour plot of \(f> + \ for the same parameter. The locations of vortices for 
(f>+ can be calculated from the phase distribution 9+(x,y) = sin _1 (Im(/)+(x, y)/\<j)+ (x, y)\). There exist a vortex at a 
point, if the integral of the phase grandient along an anticlockwise path encircling the point is a nontrivial multiple of 
2tt. We have counted the path integral by discretizing the (x, y) space with Ax = 1/64. Figure 2(c) shows positions 
of vortices of vorticity ±1 with square and x marks. The vortex cores locate near (0, —0.28) and (0, —0.48) for <fi + . 
In generic cases, there is a vortex of vorticity 1 or -1 at a position satisfying \<f>+\ — 0, where a line of Re</> + = 
intersects with a line of Im0 + = We do not show explicitly the positions of vortices later in Fig. 3 and Fig. 4, 
however, we have checked the existence of vortices of vorticity 1 or -1 at positions satisfying \<f>±\ = by calculating 
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FIG. 2: (a) |</>+| (solid curve) and \<f>-\ (dashed curve) along the line x — for A = 3n/2, and k x — 3n/2. (b) Contour plot 
of \<j>+\- ( c ) Square shows a vortex with vorticity 1 and x shows a vortex with vorticity -1. (d) Minimum values of \tf>+\ as a 
function of A for k x = A and k y = 0. 



the phase distribution. Figure 2(d) shows the minimum value of \4>+ \ as a function of A for k x — A, k y — at e = 5. 
The minimum value becomes zero and a vortex-antivortex pair appears for A > A c ~ 4.2. 

Because <f)± are periodic functions with wavelength 1, <f>± can be expressed as the simplest approximation: 

<P+ = Go++Oi+e + 2 +e + C 3+ e + C-4+e a , 

0- = Co-+Ci-e +G2-e + L 3 -e a + G4-e y . (5) 

Substitution of this ansatz into Eq. (3) yields 

mCo± = (k 2 x + k 2 y )C 0± /2-(e/2)(C 1± + C 2 ± + C 3± + C 4± ) + X(±ik x + k y )C 0T , 
fiC 1± = {(k x + 2ir) 2 + k 2 y }C 1± /2-(e/2)C ± + X{±i(k x + 27r) + k y }C lTl 
liC 2 ± = {(k x ~2Tr) 2 + kl}C2±/2~(e/2)Co± + \{±i(k x + 2n) + k y }C2 T , 
M<? 3 ± = {k 2 x + (k y + 2Tr) 2 }C 3± /2~(e/2)C ± + X{±ik x + {k y + 27:)}C 3T 

liC i± = {k 2 x + (k y -27r) 2 }C i± /2-(e/2)C ± + \{±ik x + (k y ~2TT)}C iT . (6) 
For k y = 0, Co_ = iCo+, Ci_ = iC\ +1 C 2 - — 1C2+ are satisfied, and then 
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where /1 is given by a solution of the equation 
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2 /x - (k x + 2tt) 2 /2 + X(k x + 2ir) fi - (k x - 2n) 2 /2 + X(k x - 2tt) 

e 2 2/i-(fcg+4 7 r 2 )-2Afc :r 
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(8) 



Furthermore, C3„ = iC 3 , ,64- = «C| + are satisfied. Here, * denotes the complex conjugate. The dashed curve in 
Fig. 1(a) denotes fi(k x ) by Eq. (8) at A = 7t/2. The approximation is good for A = it/2 but is not so good for large A, 
because the higher harmonics is necessary for the expansion in Eq. (5). We can assume that Co+,Ci+ and C2+ arc 
real numbers and C4+ = C| + = ReC3 + — TmC3 + . Then, 4> + and 4>- are expressed as 

4>+ = C 0+ + (Ci+ + C 2+ ) cos(2ttx) + i(C 1+ - C 2+ ) sin(27ra;) + 2ReC 3+ cos(2?n/) - 2ImC 3+ sin(2?n/), 

(p- = iC + + i(C 1+ + C 2+ ) cos(2ttx) - (C 1+ - C 2+ ) sin(27ra) + 2iImC 3+ sin(27ry) + 2iReC 3+ cos(2tt?/). (9) 
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FIG. 3: (a) Contour plot of \tp+\ for A = ir,g = 1,7 = 0.5 and L — 8. k x = k y axe evaluated as 37r/4. (b) Contour plot of 

(c) Contour plot of \<f>+\ to the linear equation Eq. (3) for A = n, k x = k y = 37r/4. (d) Minimum values of \4>+\ as a function of 

A for k x — k y — A/V2- 
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FIG. 4: (a) Contour plot of for A = n,g = 1,7 = 2 and L = 8. k x = k y axe evaluated as 37r/4. (b) Contour plot of ji/'-l- 
(c) Contour plot of the superposition of | (</>++ + 0+-)/V2| to the linear equation Eq. (3) for X — tt and k x — k y — ±37r/4. (d) 
Contour plot of the superposition of |(0++ + 0+_)/v / 2| to the linear equation Eq. (3) for A = 1 and k x — k y — ±7r/4 



Im0+ = and Re0_ =0 arc satisfied on the line x — 0. When A is small, the minimum values of Re</>+ and lm.4>- 
are positive and there are no vortices. When A is increased the minimum values decrease and reach 0, and then a 
vortex pair is created. A vortex core of <t>+ is located at a point on the line x = where Re </>+(0, y) = is satisfied, 
and similarly a vortex core of 0_ is located at a point on the line x = where Im <p_ (0, y) = is satisfied. 

Even for g and 7 is not zero, the Bloch state is a good approximation for the stationary state for 7 < g. We 
have performed numerical simulation of Eq. (1) by the imaginary time evolution method similar to Eq. (4) and found 

stationary solutions. The system size is L x x L y = L x L and the total norm N = f Q J Q (IV'+I 2 + \ip-\ 2 )dxdy is set 
to be L 2 in this paper. Periodic boundary conditions are imposed. The potential is shifted as U = — e[cos{27r(x — 
1/2)} + cos{27r(y — 1/2)}] by (1/2, 1/2) to confine the wave pattern in the range of [0, L\ x [0, L\. 

Figure 3(a) and (b) show contour plots of an d \^-\ at g = 1,7 = 0.5, L = 8, A = n, and e = 5. The contour 
plot is drawn in the region [0,2] x [0,2], and the contour lines are drawn for \ip±\ = 0.025,0.05,0.075,0.1, 1 and 1.5. 
Vortex pairs exist in each cell of size 1 for this parameter, and a vortex lattice is constructed as a whole. Vortex 
lattices were experimentally found first in rotating BECs [12| and recently in BECs under synthetic magnetic fields 
by atom- light interaction [13j . In our model equation, vortices are spontaneously created by the spin-orbit coupling. 
The wavevector (k x , k y ) is evaluated at (37r/4, 37r/4). Positions of vortex cores for ip + and ip- are mutually deviated. 
Figure 3(c) shows a contour plot of \4>+\ for the linear equation corresponding to g = 0, 7 = for k x — k y — 3ir/4 at 
A = 7r and e = 5. The eigenvalue /i takes a minimum at (k x ,k y ) = (37r/4, 37r/4) in the finite size system of L — 8, 
where k x (k y ) takes a discrete value 2im x /L (2irn y /L) with integer n x (n y ). The contour plot is almost the same as 
Fig. 3(a). It means that the Bloch wave is a good approximation for the solution to the GP equation. Figure 3(d) 
shows the minimum values of \4>+\ for the linear equation as a function of A for k x = k y = A/\/2 at e = 5. The 
minimum value becomes zero and vortices appear for A > 2.9. It is related to the existence of vortices at A = tt. 

Stripe wave states are expected to appear for 7 > g. The superposition of Bloch waves of {k xi k y ) and (— k x , —k v ) 
is a simple approximation for 7 > g. Figure 4(a) and (b) show contour plots of |-0 + | and at g = 1,7 = 2, L = 8, 
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FIG. 5: (a) Spin configuration of (s x (i,j), s y (i,j)) at A = 7T,g = 1,7 = 0.5 and L — 8. (b) Spin configuration of (s x (i,j), s y (i,j)) 
at A = 7r, g = 1, 7 = 2 and L = 8. (c) Spin configuration of s z (i,j) at A = 7r, g = 1, 7 = 2 and L = 8. 



and A = tt. The wavevector is evaluated as (k x ,k y ) = (37r/4, 37r/4) in this case, too. The contour lines are drawn 
for \ijj± \ = 0.025, 0.05, 0.075, 0.1, 1 and 1.5. Vortex cores exist in dark pointed regions. The vortex lattice structure is 
rather complicated. The circular contour lines correspond to peak regions of \ip±\. The peak regions stand in a line 
in the direction of angle — 7r/4 and the peak lines for ip + and tp_ alternates in the diagonal direction of angle 7r/4. 

Figure 4(c) shows a contour plot of a linear combination | (4>++ + <fr^ )/V%\ of two Bloch waves </>++ and 4>-\ with 

(k x ,k y ) = (37r/4, 37r/4) and (— 37r/4, — 37r/4) for the + component at A = tt. The superposition of the Bloch waves 
is a good approximation for the stationary solution to the GP equation. The superposition of two Bloch waves with 
opposite wavevectors generates a standing wave. For plane waves, the amplitude becomes zero at the nodal lines. 
The nodal lines are perturbed by the optical lattice and vortices are generated. A vortex lattice structure therefore 
appears even for small A in case of 7 > g. Figure 4(d) shows a vortex lattice pattern with k x = k y — 7r/4 at A = 1 
and e = 5. For large A, a vortex pair is created in a single Bloch wave and the superposition of the two Bloch waves 
make the vortex lattice structure more complicated as shown in Fig. 4(c). 

The complicated patterns might be simplified, if a spin representation is used, which was discussed in the Bose- 
Hubbard model [y, 0|. The whole system is divided into cell regions of [i — X [j — 1, j]. The spin variables 
s x (i, j), s y (i, j) and s z (i,j) are defined for each cell labeled by as 

s x (i,j) = j I ijy a x ipdxdy = / / (tp+ip- + iplip+)dxdy, 

s v(hj) — [ f ^ Uyipdxdy = ( ( (-iip+ ip- + ip+)dxdy, 

Ji-lJj-l Ji-lJj-l 



s z (i,j) = I [' ^a z i]>dxdy= [ [ (|V + | 2 -|^-| 2 )^, (10) 

Ji-lJj-l Ji-lJj-l 

where a x , a y and a z are the Pauli matrix, and t denotes the complex conjugate transpose. Figure 5(a) shows (s x , s y ) 
corresponding to the pattern in Figs. 3(a) and (b) for g = 1, 7 = 0.5, A = tt and e = 5. The vector (s x (i,j), s y (i, j)) is 
expressed as an arrow on each lattice point at (i — 1/2, j — 1/2). The pattern is interpreted as a ferromagnetic state 
in the (x,y) plane in this spin representation. The spin s z is zero for this pattern. Figures 5(b) and (c) show spin 
configurations respectively for (s x , s y ) and s z for the pattern at g = 1 and 7 = 2 shown in Figs. 4(a) and (b). The spin 
configuration is also rather complicated. The wavelength of the spin configuration is 4 both in the i and j directions. 
An anti- ferromagnetic order is seen in the diagonal direction of angle ir/4 and a ferromagnetic order appears in its 
orthogonal direction of angle — 7r/4 both for the (s x , s y ) and s z patterns. The (s x , s y ) component appears at the sites 
where the s z component vanishes, and the s z component appears at the sites where the (s x , s y ) component vanishes. 

To summarize, we have studied the Gross-Pitaevskii equation with spin-orbit coupling in an optical lattice. We have 
found that a vortex lattice structure appears for large A in case of 7 < g. A vortex lattice structure appears even for 
small A in case of 7 > g, because the nodal lines in the stripe wave pattern are perturbed by the optical lattice. We have 
found a complicated spin configuration in a case of 7 > g. The complicated patterns can be qualitatively understood 
by the corresponding Bloch waves. The Bloch waves are further approximated by a Fourier series expansion with five 
modes to understand the formation of the vortices. We have obtained various spin configurations by changing the 
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parameter A. The detailed phase diagrams by changing various parameters are under study. 
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